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ABSTRACT 


This  report  describes  a  program  for  simulating  the  motion  of  a  rigid 
body  containing  a  reaction  wheel,  damping  media  and  torquing  (magnetic  or 
otherwise).  The  dynamical  equations  are  simulated  in  terms  of  the  angular 
velocities  in  body-fixed  coordinates  and  the  kinematic  equations  are  simu¬ 
lated  in  terms  of  direction  cosines.  The  entire  system  of  equations  has  been 
programmed  using  a  fourth  order  predictor- corrector  method  with  automatic 
step- size  adjustment  and  a  fourth  order  Runge-Kutta  starting  routine.  The 
ratio  of  computer  time  to  real  time  is  between  1  to  lOand  2  to  1  for  LES-5- 
type  satellites.  The  potential  user  need  only  be  familiar  enough  with  this 
report  to  fill  out  the  input  data  sheets  on  pages  vii  -  ix. 


Accepted  for  the  Air  Force 

Franklin  C,  Hudson 

Chief,  Lincoln  Laboratory  Office 


CONTENTS 


Abstract  iii 

Input  Data  Sheet  vii 

I.  Introduction  1 

II,  Dynamical  Models  And  Their  Equations  of  Motion  2 

III,  Kinematics  6 

IV,  External  Torques  10 

V,  Numerical  Procedures  and  Results  10 

VI,  Sample  Problem  1:  Torque-Free  Motion  12 

VII,  Sample  Problem  2:  A  Magnetic  Attitude  Control  Law  17 

VIII,  Sample  Problem  3:  Reaction  Wheel  Spin  Down  23 

Appendix  1:  Usage  of  the  Attitude  Stabilization  Program  28 

References  40 


V 


INPUT  DATA  SHEET 


I.  Kinematic  Specifications 
Fill  in  either  A  or  B 

A.  Initial  values  of  direction  cosines: 


A(0)  = 


B.  Initial  Euler  angles: 

9  =  rad,  cp  =  rad  ^  =  rad 

Q  -  o  -  o  — - 

IL  Dynamic  Specifications 

The  i^^  moment  of  inertia  and  the  i^^  component  of  angular  velocity  are 
assumed  to  be  taken  about  the  i^^  principal  axis  of  the  satellite.  If  there  is 

a  reaction  wheel,  the  axes  must  be  numbered  so  as  to  place  the  wheel  on  the 
o  rd 

3  axis. 


Satellite  Characteristics 
Initial  value  of  angular  velocity: 

rad 

uu(0)  =  ' 


(- 


sec , 


Moments  of  inertia: 


kg  m 


’^2  = 


rad 
sec , 


,  2 
kg  m  , 


I3  = 


rad 

sec 


) 


m 


B.  Rotor  Characteristics 


.th 


If  there  is  a  spherical  damping  rotor  on  the  i  axis,  fill  in  )(.(0),  its 
initial  angular  velocity;  R.,  its  moment  of  inertia;  and  k.,  its  damping  constant. 
If  no  rotor,  write  NONE  in  the  corresponding  blanks. 


X{0) 

Ri 


(■ 


rad 

sec. 


rad 

sec. 


rad 

sec 


) 


kg  m  ,  R2=_ 

,  2 
kg  m 


kg  m 


sec 


^^2=- 


kg  m 
sec 


^3  = 


^^3  " 


kg  m 


m 


sec 


Vll 


C.  Reaction  Wheel  Characteristics 

If  the  satellite  has  a  reaction  wheel,  fill  in  Xo>  its  initial  angular 
velocity,  and  W,  ,  ,  W^,  its  moments  of  inertia.  Also  supply  a  control  law 

determining  (A  »  X  *  the  torque  of  the  reaction  wheel  on  the  satellite. 

If  no  wheel,  write  NONET. 


X 


rad 

sec 


W  =  W 

Wi 


w. 


kg  m 


III.  External  Torque 

Supply  a  control  law  determining  n,  the  external  torque  on  the  satellite, 
in  body-fixed  coordinates. 


1 


n . 


nt  m 


nt  m 


nt  m 


IV.  Termination  Criterion 

Fill  in  A  or  B,  or  both, 
nate  the  simulation. 


If  both,  the  first  criterion  to  be  met  will  termi- 


A.  Terminate  after  integrating  over  a  period  of 


sec . 


B.  Run  the  simulation  until  the  solutions  satisfy  the  following 
condition: 


V.  Title  for  Output 

Give  a  brief  title  for  the  output,  128  characters  or  less;  a  character  is 
a  letter,  a  numeral,  a  punctuation  mark,  or  a  space. 
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NUMERICAL  SIMULATION  OF  THE  ATTITUDE  EQUATIONS  OF  A  SATELLITE 
CONTAINING  DAMPING  ROTORS  AND  A  REACTION  WHEEL 


L  INTRODUCTION 

Our  objective  has  been  to  develop  a  computer  program  for  solving  the 
equations  which  describe  the  motion  of  a  satellite  about  its  center  of  mass 
and  to  get  a  feeling  for  how  accurate  the  simulation  might  be.  The  major 
decisions  which  needed  to  be  made  were: 

(i)  What  is  the  simplest  model  of  the  satellite  which  can  account  for  the 
main  features  of  the  motion.  In  particular,  how  can  we  include 
damping  of  wobble  without  introducing  too  much  complexity? 

(ii)  Which  coordinates  should  be  chosen  in  carrying  out  the  simulation-- 
Euler  angles,  quaternions,  Cayley- Klein  parameters,  direction  co¬ 
sines  or  some  other  set. 

(iii)  What  method  of  integration  should  be  chosen,  e.  g.  Runge-Kutta, 

Adams- Bashforth,  etc.  ,  and  what  step  sizes  and  error  bounds  (in 
automatic  step  size  control)  should  be  used. 

Naturally  the  decision  in  each  case  was  made  in  terms  of  what  we  felt 
would  maximize  ultimate  usefulness  of  the  program.  The  uses  we  had  in 
mind  are  exemplified  by  the  following,  all  of  which  have  been  suggested  to 
us  by  other  people  in  the  Laboratory  working  on  LES-5,-6  and-7. 

(i)  Suppose  a  rigid  body  has  within  it  a  motor  driven  reaction  wheel. 

The  power  fails  and  the  wheel  spins  down  because  of  friction.  How 
does  this  affect  the  motion  of  the  main  body  over  a  period  of  time 
long  enough  to  let  the  wheel  lose  90%  of  its  relative  angular  velocity? 

(ii)  Suppose  a  rigid  body  has  magnetic  rods  and  a  control  law  for  orient¬ 
ing  it  along  a  magnetic  field.  How  large  can  the  magnetic  torques  be 
without  causing  instability  and  how  fast  can  such  a  system  be  made 
to  re  spond  ? 

(iii)  Given  a  particular  configuration  of  gas  jets  and  a  particular  control 
law,  how  much  fuel  will  it  take  to  orient  a  satellite  from  a  given  ini¬ 
tial  state  ? 

This  report  does  not  contain  any  true  applications  in  the  sense  of  detailed 
parametric  studies  of  a  particular  problem.  Sample  problems  however  have 
been  worked  and  the  results  are  given  along  with  the  computation  times  and 
some  remarks  on  accuracy. 

From  the  point  of  view  of  the  potential  engineering  user  the  main  fact  is 
that  by  filling  out  the  input  data  sheet  (pages  vii  -  ix)  and  handing  it  to  a 
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programmer  he  can  in  a  short  time  get  a  complete  simulation  without  ever 
studying  the  internal  workings  of  the  program.  In  order  to  assure  useability 
independent  of  programming  personnel  we  have  included  a  reasonably  com¬ 
prehensive  guide  to  the  program  in  an  appendix. 

II.  DYNAMICAL  MODELS  AND  THEIR  EQUATIONS  OF  MOTION 

The  actual  mechanisms  which  account  for  wobble  damping  in  most  cases 
are  far  too  complex  to  simulate  in  detail.  Perhaps  the  simplest  model  which 
bears  resemblance  to  the  physical  situation  is  that  of  a  small  rotor  in  an 
otherwise  rigid  body  with  the  axis  of  the  rotor  along  one  of  the  principle  axes. 
If  the  center  of  mass  of  the  rotor  coincides  with  the  center  of  mass  of  main 
body  and  if  the  rotor  is  dynamically  spherical  then  further  simplifications  oc¬ 
cur,  For  the  sake  of  symmetry  one  may  also  want  to  consider  additional  ro¬ 
tors  along  the  other  principal  axes. 

A.  Case  1:  3  Spherical  Damping  Rotors 

Lets  now  consider  the  case  of  a  rigid  body  with  moments  of  inertia  1^, 
and  about  its  principal  axes.  We  further  suppose  that  the  body  contains 
sphe  rical  rotors  having  as  their  moment  of  inertia  about  any  axis  passing 

through  their  center  of  mass,  R.,  R^  and  R  respectively.  These  rotors  are 

A  ^  ^  th 

mounted  in  such  a  way  as  to  permit  the  i  rotor  to  rotate  about  the  i  prin¬ 
cipal  axis.  If  we  let  uu ,  ,  and  denote  the  angular  velocities  of  the  main 

1-  ^  j  t  h 

body  and  if  denotes  the  angular  velocity  of  the  i  rotor  (see  Fig.  1),  then 
the  full  equations  of  motion  are: 


(R^  +  R3  +  uu  j 

=  (I2- 

13)0.20)3 

+  f  ^  +  n 

(Rj+R3+l2)uu2 

^  ^3  " 

'i>  ”1  “3 

(Rj  +  1^2  +  ^3^^"  3 

=  ('1  - 

'2*  “1  *2 

Ri  Xi 

^2  2 

^  2 

R3  X  3 

■  3 

2 


l3-63-86?6| 


Fig.  1.  A  rigid  body  with  three  pivoted  spherical  rotors  at 
the  center  of  mass.  See  section  2.  1,  case  1,  for  equations 
of  motion.  Viscous  frictions  between  rotors  and  main  body 
provides  damping.  Rotors  are  shown  withdrawn  for  ease  of 
visualization. 
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where  i  is  the  component  of  the  reaction  torque  between  the  i^^  wheel  and 
the  main  body  lying  along  its  axis  of  rotation.  To  get  damping  we  can  assume 
viscous  friction  in  the  bearing  supporting  the  reaction  wheels.  If  this  is  lin¬ 
ear  then 


k.  2  0 
1 


and  the  above  equations  become 


(R2  + 

R3  + 

ij)  = 

(I2  -  I3)  ^2  ^^3  ^1  ^^1  “  '^1 

(Rj  + 

R3  + 

■2' *2  = 

(I3  —  ^3  4  "  ^2^  *^2 

(R^  + 

R2  + 

I3)  ij  = 

(I^  -  I^)  +  k3  (X3-  UU3)  +  n3 

R}  X  1 

-  X  j) 

^2^2 

•'2'"’2  -  ^2> 

R3  X  3 

^3(0)3  -  X  3) 

B.  Case  2 

A 

Reaction  Wheel  and  2  Spherical  Damping  Rotors 

We  now  consider  the  possibility  of  a  reaction  wheel  along  one  of  the 
principal  axes,  say  the  3-axis.  To  keep  the  notation  simple  we  will  suppose 
that  there  is  no  damping  rotor  about  this  axis  but  that  spherical  damping  ro¬ 
tors  may  be  present  along  the  other  two  principal  axes.  It  does  not  seem  to 

be  general  enough  to  require  the  reaction  wheel  to  be  dynamically  spherical 
so  we  let  it  be  cylindrical  with  moments  of  inertia  and  (see  Fig.  2). 

In  the  program,  the  coordinate  axes  are  always  in  such  a  position  that  ~  ^2 ' 
In  the  interest  of  simplicity  we  assume  that  the  center  of  mass  of  the  reaction 
wheel  coincides  with  the  center  of  mass  of  the  main  body.  In  this  case  the 
equations  are 

(ii  +  R2  “  ^^2 "  ^3^  ^2  ^3  ^^2 '  ^3^  ^2  ^3  -^1 

(I2  +  Rj  +  W3)  (I)  2  =  (I3  -  I^)  11)3  +  (  W3  -  W^)  uUj  X3+^  2  ^  ""2 
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Fig.  2.  Figure  1  with  the  rotor  about  the  3-axis  replaced  by  a 
cylindrical  reaction  wheel.  See  section  2.2.  The  two  rotors 
are  spherical. 
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(I3  +  Rj  +  R^)  ^1^3  -  (Ij  -  ^2)  ^2  ^3  "^3 

R,  X,  =  -<i 

^2^2  “  "  *  2 

'^3^3  =  -'3 


As  before,  f .  is  that  component  of  the  reaction  torque  of  the  i  ro- 

^  th 

tor  on  the  main  body  which  lies  along  the  axis  of  rotation  of  the  i  rotor. 

Typically  we  let  f  ^  and 

^1  "  ^1  ^^1  ■  ^1^ 

^2  ~  ^2  ^'^2 "  ^2^ 


The  choice  of  like  the  choice  of  n^,  n^  and  n^  must  be  made  so 
as  to  achieve  the  desired  attitude. 

III.  KINEMATICS 

One  has  several  choices  when  it  comes  to  describing  the  orientation  of  a 
rigid  body  with  respect  to  an  inertially  fixed  reference  frame.  As  is  well 
known,  the  popular  Euler  angles  are  not  well  suited  for  numerical  work  be¬ 
cause  of  singular  points.  The  quaternions  are  superior  in  this  respect  but 
because  of  our  need  for  direction  cosine  information  in  computing  magnetic 
torques  we  have  chosen  to  use  direction  cosines  throughout.  This  necessi¬ 
tates  integrating  9  kinematic  equations  instead  of  4,  as  would  be  the  case  for 
the  quaternions,  but  these  equations  are  well  behaved  and  for  problems  which 
would  require  the  calculation  of  direction  cosines  at  any  stage  this  seems  like 
an  appropriate  choice.  (We  have  thus  far  not  compared  integration  times 
with  an  integration  procedure  based  on  quaternions.  ) 

Figure  3  shows  a  moving  triad  and  a  fixed  triad  having  a  common  origin. 

t  h 

We  define  a.,  as  the  orthogonal  projection  of  a  unit  vector  along  the  i  axis 
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|3-<3‘8e28| 


Fig.  3.  The  moving  triad  and  the  fixed  (inertial)  triad. 
Note  definition  of  ^^2*  defined  analogously. 
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of  the  moving  system  on  the  j  axis  of  the  fixed  system.  If  we  let  OJ.  be  the 

th  ^ 

angular  velocity  of  the  moving  system  about  its  i  axis  then  the  equations  of 
evaluation  for  the  direction  cosines  are 


11 

II 

^21 

(U3  - 

12 

= 

^22 

UJ3  - 

13 

= 

^23 

^^3  “  ^33 

21 

= 

^31 

^1  "  ^11 

22 

^32 

^1 

“  ^12 

23 

= 

^33 

-  aj3 

31 

= 

^11 

^2 

“  ^21 

32 

= 

^12 

^2 

“  ^22 

33 

= 

^13 

^"2 

“  ^23 

UL) 

UL) 

0) 

OJ 


2 

2 

2 

3 


uu 

0) 

uu 

CD 

OJ 


3 

3 

1 

1 

1 


Expressed  more  succinctly, 

A  =  A 

where 


^11 

^12 

^13 

0 

^3 

-  ^2' 

^21 

^22 

^23 

3 

= 

-0)3 

0 

^1 

-  ^31 

^32 

^33  - 

- 

-0^1 

0 

We  remark  in  passing  that  A  is  an  orthogonal  matrix  so  A’A  is  the  3  by  3  iden¬ 
tity  matrix,  where  prime  denotes  transpose. 

For  the  sake  of  completeness  we  include  the  formula  relating  the  direc¬ 
tion  cosines  and  the  Euler  angles. 
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A.  Euler  Angles  in  Terms  of  Angular  Velocities  and  Direction  Cosines 


0 

cp 

6 

cp 

'1' 


arccos  a 


33 


arccos 

arccos  (a^^/sitiQ) 


UOj  cos  \|;  -  0)^  sin  i|/ 

UO^  sin  \|f/sin  0  +  \lf/sin  0 


'3  "^2  ^  9/sin  0  -  uo^  sin  ij/  cos  0/sin  0 


0)^  -  uo. 


B.  Directions  Cosines  in  Terms  of  Euler  Angles 


c  ^ 

c  cp  - 

C  0 

s  cp 

s  if 

c  if 

S  cp  +  c 

0 

c  cp 

s  'll 

s  'if 

s  6 

s  ij/ 

c  cp  - 

C  0 

S  cp 

C  'il 

-  S 

S  cp  +  C 

0 

c  cp 

C  l)j 

c  ^ 

s  0 

s  0 

S  cp 

-  s  0 

c  cp 

c  0 

C.  Angular  Velocities  in  Terms  of  Euler  Angles 


uo  j  =  cp  sin  0  sin  \|f  +  0  cos  'if 


U) 


2  =  ^  sin  0COS  -  0  sin 


UO^  =  cp  cos  0  +  ij/ 

In  many  calculations  the  angular  momentum  plays  a  key  role.  It  may 
be  shown  that  the  matrix 

*  -  1  -  1 

L  =  (A^  _Q^  A)det^ 

is  antisymmetric  and  -^23’  ”^13  ^12  components  of  the  angular  mom¬ 

entum  along  the  1*,  2*  and  3*  axes,  respectively,  in  the  inertial  triad. 


is  of  some  interest  to  note  that  the  Euler  equations  for  a  rigid  body  can  be 
written  in  this  matrix  notation  as 


(det  I)  l"^  D  l"  ^  I  -  I  +  N 
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IV.  EXTERNAL  TORQUES 


It  is  essential  to  keep  in  mind  that  the  basic  dynamical  equations  are  writ¬ 
ten  in  body  fixed  coordinates.  The  torques  must  be  expressed  in  body  fixed  co¬ 
ordinates  and  must  be  resolved  along  the  principal  axes  of  the  body.  Thus  our 
equations  take  the  simplest  form  if  thrusters,  magnetic  rods,  etc.  are  located 
in  such  a  way  as  to  make  the  resolution  of  torques  along  the  principal  axes  sim¬ 
ple.  Because  of  the  importance  of  magnetic  torquing  we  will  illustrate  the  cir¬ 
culation  of  torques  using  this  as  an  example. 

Consider  a  magnetic  field  of  intensity  H  directed  along  the  3-axis  of  a  co¬ 
ordinate  system  fixed  in  space  .  Suppose  there  is  a  set  of  3  orthogonal  magnets 

having  moments  m,,  m^  and  m^,  respectively.  To  calculate  the  torque  about 
th  ^  ^  th 

the  i  moving  axis  due  to  the  j  magnet  we  needto  compute  the  cross  product  of 
the  two  vectors.  Expressed  in  terms  of  the  moving  coordinate  system,  h  has 
the  components 


^11  ^12  ^13 

"o  ■ 

"^3‘ 

^21  ^22  ^23 

0 

=  h 

^23 

-^31  ^32  ^33- 

_  h  , 

-^33- 

Hence  we  have 


■^r 

-mi" 

^13 

"^2 

^33 

-  ^3 

^23' 

^2 

^2 

X 

^23 

m3 

^13 

- 

^33 

-'^3- 

-”^3- 

-^33- 

-^1 

^23 

^2 

^13- 

V.  NUMERICAL  PROCEDURES  AND  RESULTS 

The  attitude  stabilization  program  employs  a  fourth-order  predictor- 
corrector  method  of  numerical  integration,  with  the  predictor  calculated  from 
the  Adams- Bashforth  formula,  and  the  corrector  from  the  Adams- Moulton  for¬ 
mula.  Because  iterating  the  corrector  formula  proved  very  costly  in  execution 
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time  and  only  marginally  profitable  in  error  reduction,  the  first  approximation 
to  the  corrector  formula  has  been  taken  as  the  final  value.  Starting  values  for 
the  Adams-Moulton  routine  are  obtained  by  fourth-order  Runge-Kutta  integra¬ 
tion. 

To  achieve  the  greatest  accuracy  in  the  least  time,  careful  control  must 
be  exercised  over  the  integration  step  size.  Therefore,  after  each  step,  a 
measure  of  the  truncation  error  is  evaluated  as  follows  for  each  of  the  fifteen 
integral  curves. 

1  =  I  (Ve  -  yp)/ye  I 

where  y  is  the  corrected  value  of  y,  and  y^  is  the  predicted  value.  (Since  the 

quotient  may  be  misleading  when  y  is  near  zero,  T\  is  set  to  zero  whenever 

jy  I  ^  5  X  10  .  )  The  step  size  for  the  next  step  is  determined  by  the  value 

of  T|  ,  the  largest  of  the  fifteen  error  estimates.  Lower  and  upper  bounds 
max  ^ 

on  Tj  ,  GLB  and  LUB,  are  supplied  as  input  data.  If  T]  is  less  than 
max  max 

GLB,  the  step  size,  h,  is  doubled;  if  Tl  is  greater  than  LUB,  h  is  halved. 
When  the  step  size  is  changed,  the  new  starting  values  are  computed  by  fourth- 
order  Runge-Kutta  integration. 

Because  the  Runge-Kutta  method  requires  twice  as  many  derivative  evalu¬ 
ations  per  step,  it  was  expected  to  be  about  twice  as  slow  as  A.dams- Moulton 
integration.  As  it  turned  out,  an  Adams- Moulton  step  without  step  size  control 
was  only  1.2  times  as  fast  as  a  Runge-Kutta  step,  and  an  Adams- Moulton  step 
with  step  size  control  took  about  as  much  computer  time  as  Runge-Kutta  step 
did.  Although  nearly  as  fast  in  our  application  as  the  Adams- Moulton  method, 
Runge-Kutta  integration  is  not  amenable  to  step  size  control.  Therefore, 
Runge-Kutta  has  been  relegated  to  the  starting  routine,  and  Adam s- Moulton 
integration  has  been  chosen  for  the  actual  simulation. 

Double  precision  (about  16  significant  decimal  digits)  is  used  throughout 
the  program,  for  single  precision  (about  8  significant  digits)  proved  to  be  in¬ 
adequate.  After  only  30  sec  of  simulation  time,  the  discrepancy  between  the 
UD  integral  curves  obtained  with  single  vs.  double  precision  is  already  0.2% 
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even  with  a  step  size  as  small  as  0.01  sec.  The  most  commonly  observed 
mean  step  sizes  are  about  ten  times  as  big,  and  as  step  size  increases,  the 
gap  between  single  and  double  precision  accuracy  widens  exponentially. 

-7  -  5 

With  GLB  =  10  and  LUB  =  8  x  10  ,  the  program  has  the  following  speed 

and  accuracy.  The  ratio  of  simulation  time  to  real  time  is  on  the  order  of  10:1, 
with  the  exact  value  determined  by  the  complexity  of  the  derivative  functions. 
Over  the  first  two  hours  of  simulation  time,  the  mean  error  in  the  solutions 
does  not  exceed  0.2%;  over  the  first  three  hours  it  does  not  exceed  0.  3%.  The 
error  may  remain  well  below  these  upper  limits. 

VI.  SAMPLE  PROBLEM  1:  TORQUE- FREE  MOTION 

In  order  to  fix  ideas  and  to  check  the  numerical  procedure  against  a  solu¬ 
ble  problem,  we  consider  first  the  case  of  the  torque-free  motion  of  a  spinning 
cylindrical  body  without  damping.  The  physical  situation  is  shown  in  Fig.  4. 

The  completed  data  sheets  are  on  the  following  two  pages.  In  order  to  get  an 
analytical  solution  we  observe  that  the  dynamical  equations  are 

=  [(I  -  ^3 

^2  =  C  (I3  -  I  )/l]  ^3 

0)3  =  0 

Thus  letting  0)  =  03^(0)  (I  -  l^)/l  we  have  as  a  solution 


p  -I 

.  — 

—  — 

% 

cos  UU  t  sin  0)  t  0 

u)^(O) 

-  sin  ua  t  cos  o)  t  0 

0^(0) 

1 - 

0  0  1 

_  11)3(0)  _ 

If  we  select  the  initial  time  in  a  suitable  way  then  ®  kinematical 

equations  are  A  =  Q  A  with  given  by 
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INPUT  DATA  SHEET 


1.  Kinematic  Specifications 
Fill  in  either  A  or  B 

A.  Initial  Values  of  Direction  Cosines 


A(0) 


73/2 


B.  Initial  Euler  Angles 


0  = 
o 


rad,  cp  = 
-  o 


rad, 


rad 


II.  Dynamic  Specifications 

The  i^^  moment  of  inertia  and  the  i^^  component  of  angular  velocity  are 
assumed  to  be  taken  about  the  i^^  principal  axis  of  the  satellite.  If  there  is 

a  reaction  wheel,  the  axes  must  be  numbered  so  as  to  place  the  wheel  on  the 
3rd 


axis. 


A.  Satellite  Characteristics 

Initial  value  of  angular  velocity: 

rad 
0  sec. 


ui)(0)  = 


(■ 


6/3  rad 


Moments  of  inertia: 


15 


sec, 


rad 

sec 


) 


5  kg  m  ,  I2  = 


5  kg  m  ,  I,  = 


B.  Rotor  Characteristics 


th 

If  there  is  a  spherical  damping  rotor  on  the  i  axis,  fill  in  X-(0)>  its 
initial  angular  velocity;  Rj,  its  moment  of  inertia;  and  k.,  its  damping  con¬ 
stant.  If  no  rotor,  write  NONE  in  the  corresponding  blanks. 

/  rad  rad  rad  ^ 


R, 


=  I  NONE 

sec, 

NONE 

sec 

\ 

=  NONE 

,  2 

kgm  , 

^2  = 

NONE 

,  2 

kgm  , 

=  NONE 

,  2 

sec  , 

=  , 

NONE 

sec  , 

NONE 


sec 


k3  = 


2  ' 

NONE 

kgm 

u  2 
kgm 

NONE 

sec 
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C.  Reaction  Wheel  Characteristics 


If  the  satellite  has  a  reaction  wheel,  fill  in  X  3*  initial  angular 
velocity,  and  Wj,  W2,  W3,  its  moments  of  inertia.  Also  supply  a  control  la 
determining  l^(A,  X3)j  the  torque  of  the  reaction  wheel  on  the  satellite, 
no  wheel,  write  NONE. 

rad 

X3  =  NONE  sec 
=  NONE  kgm^, 


W, 


NONE  kg  m 


i 3  =  NONE 


III.  External  Torque 

Supply  a  control  law  determining  n,  the  external  torque  on  the  satellite, 
in  body-fixed  coordinates. 


II 

0 

nt 

m 

II 

0 

nt 

m 

II 

0 

nt 

m 

IV.  Termination  Criterion 

Fill  in  A  or  B,  or  both.  If  both,  the  first  criterion  to  be  met  will  termi¬ 
nate  the  simulation. 

A.  Terminate  after  integrating  over  a  period  of  10,  800 _  sec. 

B.  Run  the  simulation  until  the  solutions  satisfy  the  following  condition: 
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l3-63-n29l 


Fig.  4.  Sample  problem  Number  1.  Torque  free 
motion 
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n(t)  = 


-0)3(0) 


-  0)^(0)  sin  uut 
0)  (0)  cos  (A)t 


0)^(0)  sin  UJt  -W^(0)  cos  tot 


as  can  be  readily  inferred  from  the  classical  treatments 


COS  (JDt 

sin 

cut 

0 

1 

0 

0 

cos 

cp  t 

0 

sin  cp  t 

0 

0 

A(t)  = 

-  sin  UUt 

cos 

cut 

0 

0  cos 

0  sin 

0 

-  sin 

cp  t 

0 

cos  ^  t 

0 

0 

0 

0 

1 

0-  sin 

0  cos 

0 

0 

0 

0 

where 


;os  0  =  1_  UU_/ 


"2  2  2  2  2  2 

1  to^  +  1  to^  +  I3  0)3 


) 


4 


2 ,  2 
"1  +”2 


sin  0 


In  the  case  of  torque-free  motion,  the  Adams- Moulton  solution  curves 
can  be  compared  directly  with  the  exact  solutions.  Let 


p.j(t)  =  |(a.^.(t)  -  a_(t)  )/ a.j  (t) 


where  a..(t)  is  the  ij^'^  element  of  A(t)  as  computed  by  Adams- Moulton  inte¬ 
gration,  and  .  (t)  is  the  corresponding  element  as  found  by  evaluating  the 
expression  for  the  exact  solution.  Then  calculate  the  mean  value 


^a(‘> 


and  the  variance 


4  s 


1.  j 


aA(t)  = 


I 


IJ 


A(t)^ 


1.  J 
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Finally,  let  T]-  (t)  be  the  error  in  1.  -(A,t),  a  component  of  angular  momentum 

“^13  . 

computed  from  the  integral  curves  for  A.  Since  the  angular  momentum  is  con¬ 
stant  in  torque-free  motion 


(t)  =  |(i^3(A,t)  -  1^3(A,0))/i^3(A,0)  I 

2 

The  values  of  p  ,  O  ,  and  X].  tabulated  below  were  calculated  after  integrating 
A  A 

with  an  initial  step  size  of  0.  08  sec  and  with  the  truncation  error  in  each  function 

-  7  “5 

confined  to  the  interval  [lO  ,  8x10  ].  The  ratio  of  simulated  time  to  real 

time  was  10.  3  to  1. 


t 

Ka(‘) 

(A) 

13 

3600  sec 

0. 0012 

2. 0  X  10"^ 

0. 0011 

7200  sec 

0.  0021 

4.  4  X  10"^ 

0. 0022 

10800  Sec 

0.  016 

1.  0  X  10"^ 

0. 0033 

The  anomalous  value  of  p^(10,  800)  is  more  likely  an  indication  of  high  er¬ 
ror  in  the  calculated  ’’exact”  solution  than  of  high  error  in  the  integral  curves. 

2 

This  conclusion  is  supported  by  the  high  value  of  800)  and  the  low  value 

of  ri-  (10,  800).  The  preceding  values  of  p  .  may  also  be  considerably  higher 

“^13 

than  the  actual  integration  errors.  In  fine,  with  the  initial  step  size  and  error 
bounds  given  above,  and  with  an  integration  period  of  three  hours,  the  mean  error 
in  the  nine  direction  cosine  solution  functions  increases,  in  all  likelihood,  by  no 
more  than  0.  1%  per  hour. 

Since  the  expressions  for  the  exact  solutions  for  ^  do  not  contain  products 
of  trigonometric  functions  of  t,  the  error  in  these  functions  is  not  compounded. 
Moreover,  the  integrated  values  of  ^  were  in  perfect  agreement  with  the  exact 
solutions;  during  the  entire  three  hours  of  simulated  time,  there  was  no  detect¬ 
able  error  in  ^ .  Since  the  integrated  and  the  calculated  solutions  were  known  to 
the  nearest  ,  00001  sec  \  the  maximum  absolute  error  in  ^  was  5x10^  sec 
VII.  SAMPLE  PROBLEM  2:  A  MAGNETIC  ATTITUDE  CONTROL  LAW 

We  now  illustrate  how  the  program  might  be  used  to  check  out  an  attitude 
control  system.  We  assume  a  spinning  cylindrical  body  with  two  damping  ro¬ 
tors  and  no  reaction  wheel.  We  also  assume  magnetic  torquing  rods  perpendi¬ 
cular  to  each  other  and  the  spin  axis.  For  simplicity  we  examine  a  fixed 
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magnetic  field  which  is  parallel  with  the  desired  orientation  of  the  spin  axis. 
The  dynamical  equations,  a  synthesis  of  the  basic  equations  of  Sections  II,  III 
and  IV  are 


where 


(I  +  2R)  0)^ 

— 

(I  -  I3)  0)3  a)^  -  k(  0) ^ 

-  Xj)  +  n 

(I  +  2R) 

= 

(I3  -  I)  0)3  uJj  -  k(  0)2 

-  X2)  +  n 

(2R  +  13)0)3 

^3 

R  Xj 

k  (0)^  -  X;^) 

R  ^2 

k  (0^  -  X2) 

U2  a33 

*^2 

= 

1  ^33 

*^3 

^  1  ^23  "  ^2^13 

and  the  u's  (proportional  to  the  strength  of  the  respective  magnets)  are  given 
by  the  control  law 

Uq  ^^33  ^^22^2  "  ^12 

Uq  sgn  [a33  ^12 

That  is,  the  u^s  are  chosen  in  such  a  way  as  to  drive  the  angular  momentum 
vector  to  the  vertical.  The  numerical  constants  can  be  found  on  the  input  data 
sheets. 

A  plot  of  the  angle  between  the  vertical  and  the  axis  is  shown  for  two 

particular  values  of  u  - 
^  o 
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INPUT  DATA  SHEET 


I.  Kinematic  Specifications 
Fill  in  either  A  or  B 

A.  Initial  Values  of  Direction  Cosines 


A(0)  = 


0 

0 


_ 0_ 

V^/2 

-1/2 


_ 0_ 

1/2 

VT/Z 


B.  Initial  Euler  angles 


e 


rad,  cp  = 
-  o 


rad,  \lf  = 
-  o 


rad 


IL  Dynamic  Specifications 

The  i^^  moment  of  inertia  and  the  i^^  component  of  angular  velocity  are 
assumed  to  be  taken  about  the  i^^  principal  axis  of  the  satellite.  If  there  is 
a  reaction  wheel,  the  axes  must  be  numbered  so  as  to  place  the  wheel  on  the 
3^^  axis. 


Satellite  Characteristics 

Initial  Value  of  Angular  Velocity 

rad 

0  sec. 


uu  (0)  = 


(- 


rad 

sec. 


rad 

sec 


) 


Moments  of  Inertia: 


m 


^2  = 


5  kg  m  ,  1^=  6  kg 


m 


B.  Rotor  Characteristics 


.th 


If  there  is  a  spherical  damping  rotor  on  the  i  axis,  fill  in  X:{0),  its 
initial  angular  velocity;  R. ,  its  moment  of  inertia;  and  k.,  its  damping  con¬ 
stant.  If  no  rotor,  write  NONE  in  the  corresponding  blanks. 


X(0) 

^1 


( 


0 

rad 

sec, 

0 

rad 

sec. 

rad 

NONE  sec 

.  05 

,  2  „ 
kg  m  ,  Rt  = 

.  05 

,  2 
kg  m  , 

R,  =  NONE  kg 

w  ..  .  .  ^ 

1  2 
kg  m 

1  2 

^  - a 

2 

kg  m 

.  01 

sec  ,  k^  = 

.  01 

sec  , 

k^  =  NONE  sec 

) 
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C.  Reaction  Wheel  Characteristics 


If  the  satellite  has  a  reaction  wheel,  fill  in  Xo>  its  initial  angular 
velocity,  and  W^,  W^,  its  moments  of  inertia.  Also  supply  a  control 

law  determining  (A,  X the  torque  of  the  reaction  wheel  on  the  satel¬ 
lite,  If  no  wheel,  write  NONE, 


rad 

Xq 

NONE 

sec 

=  NONE  kR  =  NONE  kg 


a 


3 


NONE 


HI,  Ehcternal  Torque 


Supply  a  control  law  determining 
in  body-fixed  coordinates.  Let  Uj^  =  u 
sgn  (  (a  J  J  ^2^  ^33^ 


n,  the  external  torque  on  the  satellite, 
^  sgn  (  (a22  ^2  “  ^12  ^1^  ^33^  ^2  ~ 


= 

^2  ^33 

nt  m 

■^1  ^33 

nt  m 

— 

^1^23  ”  ^2  ^13 

nt  m 

where  u  is  either  ,01  or  ,002 
o 

IV.  Termination  Criterion 

Fill  in  A  or  B,  or  both.  If  both,  the  first  criterion  to  be  met  will  termi¬ 
nate  the  simulation, 

A.  Terminate  after  integrating  over  a  period  of  1500  sec. 

B.  Run  the  simulation  until  the  solutions  satisfy  the  following  condition: 
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Fig.  5.  Sample  problem  Number  2.  Magnetic  torquing  with  a 
fixed  magnetic  field.  Damping  rotors  are  included. 
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VIIL  SAMPLE  PROBLEM  3:  REACTION  WHEEL  SPIN  DOWN 

Consider  a  cylindrical  body  whose  moment  of  inertia  along  the  axis  of 
symmetry  is  smaller  than  the  moment  of  inertia  along  the  other  two  princi¬ 
ple  axes.  Such  a  body  might  be  held  relatively  stable  by  putting  a  spinning 
wheel  along  the  axis  of  symmetry  with  an  electric  motor  drive  to  overcome 
friction.  We  want  to  investigate  what  happens  if  the  power  fails  and  the  wheel 
spins  down  due  to  friction.  The  main  item  of  interest  might  be  the  angle  that 
the  axis  of  symmetry  makes  with  the  vertical  after  a  certain  period  of  time. 
Naturally  if  the  system  is  perfectly  balanced  and  ^  initially  then 

momentiim  lost  by  the  wheel  will  be  taken  up  by  the  body  but  no  change  in  the 
orientation  of  the  spin  axis  will  occur.  However  if  the  initial  values  of  and 
U)^  are  slightly  different  from  zero  then  wobble  will  build  up.  We  assume  the 
existence  of  damping  rotors  on  the  1-  and  2-axes  and  take  the  reaction  wheel  along 
the  3-axis.  Let  the  moments  of  inertia  be  1^  ~  ^2  ”  ^  ^3  ”  The  initial 

velocity  of  the  main  body  is  zero  and  is  assTimed  to  be  perfectly  oriented  (A(0) 
is  the  identity  matrix).  Since  the  complete  input  data  sheet  is  attached,  we 
only  remark  that  the  initial  angular  momentiim  of  the  wheel  is  180  kg  m  /sec. 

The  value  of  damping  picked  is  such  that  the  wheel  spin  down  has  a  time  con¬ 
stant  of  18  minutes.  Attached  graph  illustrates  the  misorientations  of  the  spin 
axis  as  a  function  of  time. 

When  the  reaction  wheel  spin  down  was  simulated  with  an  initial  step  size 
of  0.  04  ^^ec  and  with  the  same  error  bounds  as  in  the  torque-free  motion  c^irnu- 
lation,  the  ratio  of  simulation  time  to  real  time  was  1 1  to  1. 
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INPUT  DATA  SHEET 


I.  Kinematic  Specifications 
Fill  in  either  A  or  B 

A.  Initial  Values  of  Direction  Cosines: 


A(0) 


1 


0 


0 


0 


1 


0 


0 


0 


1 


B.  Initial  Euler  Angles: 

0 


rad,  cp  = 
-  o 


rad , 


rad 


II.  Dynamic  Specifications 

t  h.  t  h 

The  moment  of  inertia  and  the  i  component  of  angular  velocity  are 
assumed  to  be  taken  about  the  i^  principal  axis  of  the  satellite.  If  there  is  a 
reaction  wheel,  the  axes  must  be  numbered  so  as  to  place  the  wheel  on  the 
3^^  axis. 


A.  Satellite  Characteristics 

Initial  value  of  angular  velocity: 
/  rad 

UJ(0)  =  (  ,  0 5  sec,  _ . 


05 


rad 

sec, 


Moments  of  Inertia: 


5  kg  m 


5  kg  m 


rad 

sec 


) 


Jl£. 


m 


B.  Rotor  Characteristics 


.th 


If  there  is  a  spherical  damping  rotor  on  the  i  “  axis,  fill  in  X*(0)>  i-^s 
initial  angular  velocity;  R.,  its  moment  of  inertia;  and  k.,  its  damping  con¬ 
stant.  If  no  rotor,  write  ^JONE  in  the  corresponding  blaViks, 


X(0) 

^1 


rad 

sec. 


rad 

sec, 


NONE 


.  05 

.  2 
kg  m  , 

^2  = 

.  05 

,  2 
kg  m  , 

.  01 

.  2 
kg  m 

sec  , 

^2  = 

.  01 

I  2 

kg  m 

sec  , 

rad  \ 
sec  J 


sec 
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C.  Reaction  Wheel  Characteristics 


If  the  satellite  has  a  reaction  wheel,  fill  in  its  initial  angular  veloc¬ 
ity,  and  W,,  ^2*  ^3*  moments  of  inertia.  Also  supply  a  control  law  de¬ 
termining  ( A  ,  X^),  the  torque  of  the  reaction  wheel  on  the  satellite.  If  no 

wheel  write  NONE. 

rad 

X^  =  60 _ sec 

W^=W^  =  .1  kg  m^, 

^3  =  (X3  -  ^3)  X  10'^ 

III.  External  Torque 

Supply  a  control  law  determining  n,  the  external  torque  on  the  satellite, 
in  body- fixed  coordinates. 


*^3 


0 

nt 

m 

0 

nt 

m 

0 

nt 

m 

IV.  Termination  Criterion 


Fill  in  A  or  B,  or  both.  If  both,  the  first  criterion  to  be  met  will  termi¬ 
nate  the  simulation. 

A.  Terminate  after  integrating  over  a  period  of  300  sec. 

B.  Run  the  simulation  until  the  solutions  satisfy  the  following  condition: 
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Fig.  6.  Wheel  spin  down  problem.  Damping  rotors  on 
1  and  2  axes;  reaction  wheel  on  3  axis.  Initially  the  rigid 
body  has  angular  velocities  =  .05  =  .05  =  0. 
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APPENDIX  1 


USAGE  OF  THE  ATTITUDE  STABILIZATION  PROGRAM 
Purpose  of  the  Program 

The  purpose  of  the  attitude  stabilization  program  is  to  simulate  the  mo¬ 
tion  of  a  satellite  about  its  center  of  mass.  The  motion  is  described  by  a 
system  of  differential  equations  of  the  form: 


where  uu  = 


A  =  a 


.ll(t) 

=  fj(t,  A) 

.12(t) 

=  i^(t,  W,  A) 

.13(t) 

=  £^(1,  w,  A) 

.2i(t) 

=  f^(t,  A) 

=  f^(t,  UJ,  A) 

-23(t) 

=  f^(t,  W,  A) 

.3l(t) 

=  f^(t,  UJ,  A) 

.32(t) 

=  fg(t,  A) 

■33(t) 

=  fg(t,  UJ,  A) 

il(t) 

~  ^10^^’  — ’  — ’ 

A) 

i2(t) 

=  lu,  x» 

A) 

i3(t) 

A) 

:i(t) 

=  £j3(t,  X) 

.2(t) 

=  X) 

.3(t) 

~  ^1 5  (^>  X) 

Jj,  UJ2’ 

UJ3).  X  =  (Xi- 

X2. 

11 

^  12  ^13  \ 

21 

^22  ^23  I 

31 

a  a  / 

32  ^33 
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More  concisely,  in  vector  notation,  the  system  can  be  written: 


i(t) 

=  £(t,y),  where  £  =  (f  •  • 

•  ’  ^33’ 

.  .  ,  X3).  and  y  =  (a^  •  • 

•  ’  ^33’  •  •  •  >  X  3)- 

The  simulation  is  accomplished  by  using  techniques  of  numerical  integra¬ 
tion  to  compute  the  solution  curves  of  the  above  system:  A,  and  X* 

Numerical  Integration 

Given  T|,  t^,  and  the  differential  equation  y(t)  =  f(t,y),  there  is  one  and 
only  one  function  y  such  that  dy/dt  =  y  and  y(t^)  =  T].  By  numerical  integra¬ 
tion,  one  can  find  a  good  approximation  to  y,  the  solution  function,  on  a  set 
of  discrete  points,  {t^}.  The  initial  value,  T|  ,  tells  where  y  is  located  at  a 
certain  time,  t^,  and  the  value  of  the  derivative  at  this  point,  f(t^,r|),  indi¬ 
cates  where  y  is  headed.  Therefore,  if  h  is  a  small  time  interval,  one  can 
compute  with  little  error  where  y  will  be  at  t^  ~  ^  Then  y(tj^)  can  be 

used  to  evaluate  f  (t  ^ ,  y(t  ^ )) ,  which  in  turn  will  point  the  way  to  y  at  t^  =  t  +  h. 
The  process  of  getting  from  y(tj^)  to  called  an  integration  step,  and 

the  step  size  is  h  =  t.  ,  ,  -  t.. 

A  method  of  numerical  integration  is  called  self- starting  if  the  only  input 
data  required  is  t^  and  T].  The  Adam s- Moulton  predictor- corrector  method, 
which  is  used  for  the  attitude  stabilization  simulations,  is  not  self- starting. 
Only  one  value  of  y  can  be  computed  from  t^  to  T],  but  the  Adams- Moulton  for¬ 
mulas  for  y(t^_|.  require  values  of  y  at  four  consecutive  points:  t^ 
t.  .  ,t..  Therefore,  the  program  includes  a  Runge-Kutta  integration  rou- 

tine,  which  computes  the  starting  values  needed  for  the  Adams- Moulton  inte¬ 
gration. 

Runge-Kutta  Integration  Routine 

A  Runge-Kutta  integration  step  has  the  form  y(t  +  h)  =  y(t)  +  h  *  ^(t), 
where  $  (t)  is  calculated  as  follows: 
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Gl(t)  =  f(t,y)  =  y(t) 

Zl(t)  =  y(t)  +  0.  5h  Gl(t) 

G2(t)  =  f(t  +  0.  5h,  Zl) 

Z2(t)  =  y(t)  +  0.  5h  G2(t) 

G3(t)  =  f(t  +  0.  5h,  Z2) 

Z3(t)  =  y(t)  +  h  G3(t) 

G4(t)  =  f(t  +  h,  Z3) 

$(t)  =  (Gl(t)  +  2G2(t)  +  2G3(t)  +  G4(t))/6 

When  a  system  of  equations,  y(t)  =^(t,]^),  is  being  integrated,  each  func¬ 
tion  y^  in  y  has  its  corresponding  Glj,  Zly  .  .  .  ,G4yiy  which  are  computed 
exactly  as  above,  using  L  to  compute  GK,  .  .  .  important  to  note 

that  all  components  of  the  vector  Z_1  must  be  calculated  before  any  of  the  com¬ 
ponents  of  G2,  all  of  Z2  before  any  of  G3,  all  of  Z3  before  any  of  G4,  and  all 
of  before  any  of  the  next  G  1.  These  are  the  only  restrictions  on  the  order 
of  computations. 

In  the  A,  X  system,  t  appears  implicitly,  as  in  but  never  explicitly, 
as  in  sin  t,  e^,  or  3t  +  10,  Therefore,  the  Runge-Kutta  equations  reduce  to  a 
simpler  form  for  this  system; 


G1 

f(y) 

Zl 

= 

y  +  0.  5h  ❖  G1 

G2 

= 

f(Zl) 

Z2 

= 

y  +  0.  5h  G2 

G3 

- 

f(Z2) 

Z3 

y  +  h  =1'  G3 

G4 

f(Z3) 

$ 

(G1  +  2G2  +  2G3  +  G4)/6 

y 

= 

y  +  h  $ 

In  the  following  summary  of  the  Runge-Kutta  routine,  and  in  the  summar¬ 
ies  of  every  subroutine  described  below,  all  CALL’s  will  be  discussed,  and 
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all  arguments  and  COMMON  variables  will  be  listed,  together  with  all  type, 
length,  and  dimension  specifications  not  implied  by  FORTRAN  convention; 
since  no  IMPLICIT  statements  appear  in  the  program,  no  nonstandard  con¬ 
ventions  pertain.  None  of  the  COMMON  blocks  is  used  for  subroutine-to- 
subroutine  communication;  they  are  used  only  for  initializing  subroutine  vari 
ables  to  values  set  in  the  main  program,  WCHA.  Variables  labeled  ’’input” 
contain  data  supplied  to  the  subroutine  in  question;  those  labeled  ’’output”  con 
tain  data  computed  by  subroutine. 

The  Runge-Kutta  subroutine,  RK(H,  T,  Y,  FY),  computes  the  starting  val¬ 
ues  for  Adams- Moulton. 


A.  Arguments 

H,  REAL  8;  input;  step  size 

T,  REAL  8;  as  input,  T  =  t^  y  as  output,  T  =  t^. 
Y(15),REAL  8;  as  input,  Y  =  Y.(t^  output, 

Y  =  The  functions  are  ordered  as  above, 

in  the  description  of  A,  UU,  X  system. 

FY(15,  5),  REAL  8;  output;  FY( J,  K)  =  fj(t.  K)’— 
J  =  1,  15,  K  =  1,  3,  where  fj  is  the  derivative  of 
Yy  FY(J,K)  =  0,  J  =  1,  15,  K  =  4,  5. 

B.  COMMON 

/CHID0T/YD0T(15) 

REAL  *  8:  YDOT 


Adams- Moulton  Integration  Routine 


For  the  differential  equation  y  (t)  =  f(t,y),  an  Adams- Moulton  integration 
step  has  the  form: 


1 


y(t 

+ 

h)  = 

y(t)  +  h 

p(t,  f) 

f 

p 

= 

f(t  +  h, 

y(t  +  h)) 

DO 

1 

m  = 

1,  n 

y(t 

+ 

h)  = 

y(t)  +  h 

>1'  c(t,  f,  f 

f 

p 

f(t  +  h, 

y(t  +  h)) 
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where  the  functions  p  and  c  are  given  by  the  predictor  and  corrector  formulas 
described  below. 

Increasing  n,  the  number  of  iterations  of  the  corrector  formula,  has  an  ef¬ 
fect  and  a  side  effect:  the  error  decreases  but  the  computation  time  increases. 
On  account  of  the  gravity  of  the  side  effect,  n  =  1  in  the  Adams- Moulton  sub¬ 
routine;  the  first  corrected  value  is  taken  as  y(t  +  h). 

For  g,  an  arbitrary  function  of  time,  let  g^  =  g  at  t^,  n  =  i  -  3,  i  -  2,  .  .  .  , 
i  +  1.  Then  the  predictor  formula  can  be  written 

.  =  y.  +  (55  f.  -  59  f.  ,  +  37  f.  _  -  9  f.  J  *  h/24 

^1+1  1  1-1  1-2  1-3  ' 

The  predicted  value  of  y.  ,  ,  is  used  to  calculate  f .  ,  ,  ,  which  appears  in  the 
^  ^1+1  i+l 

corrector  formula: 


y 


i+  1 


y.  +  (9f 


i+  1 


+  19f.  -  5f. 

1  I 


j  +  f.  .  2)  *  h/24 


To  apply  the  Adams- Moulton  method  to  a  system  of  equations,  y  =  f  (t,  y), 

simply  apply  the  predictor  and  corrector  formulas  to  each  component  of 

making  sure  to  compute  the  entire  predicted  vector  ^  before  starting  to 

calculate  f.  ,  ,  for  the  corrector  formula. 

—1+1 

It  is  important  to  maintain  an  optimum  step  size  during  the  integration. 

If  h  is  too  small,  the  computation  takes  too  much  time,  but  if  h  is  too  big, 

there  is  a  large  error  in  the  solutions.  To  optimize  h,  an  error  estimate,  E, 

is  calculated  for  each  of  the  fifteen  functions  after  every  integration  step. 

E  ,  the  greatest  of  the  fifteen  errors,  has  upper  and  lower  bounds  speci- 
max  ^  r'r'  r- 

fied  as  input  to  the  program.  If  E  is  greater  than  LUB,  the  upper  bound, 

xnax 

then  h  will  be  halved,  and  if  E  is  less  than  GLB,  the  lower  bound,  then  h 

max 

will  be  doubled. 

Whenever  the  step  size  is  changed,  the  Runge-Kutta  starting  routine  must 
be  called,  for  the  Adams- Moulton  routine  cannot  find  the  complete  set  of  deri¬ 
vatives  at  the  new  points  {t.  i}  =  {t.-3h,  t.  -  2h,  t.-h]  unless  the 

1  1  I 

step  size  has  remained  constant. 
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The  choice  of  the  lower  and  upper  error  bounds,  GLB  and  LUB,  will  be  a 

compromise  between  accuracy  and  speed.  Often  a  good  set  of  bounds  is 
-7  -  5 

GLB  =  10  ,  LUB  =  10  .  However,  an  optimum  set  of  bounds  for  a  particu¬ 

lar  A,  X  system  can  be  determined  as  follows.  Form  n  sets  of  bounds, 

B, ,  .  .  .  ,  B  ;  B.  =  [GLB,  LUB}.  ,  i  =  1,  .  .  .  ,  n.  Run  n  simulations,  identical  ex- 
cept  in  the  choice  of  B^,  and  obtain  n  corresponding  sets  of  solutions,  , 

S^;  =  {  A,  X  i  =  Then  study  the  discrepancy  between  and  Sy 

for  {i,j  :  i  /  j  }•  Suppose  B^  has  more  restrictive  bounds  than  B^;  i.e.,  h  will 
halve  more  easily  and/or  double  less  easily  under  B^,  Then  if  the  discrepancy 
between  and  S.  is  great,  one  may  reject  Bj  as  unacceptably  lax.  If,  on  the 
other  hand,  the  discrepancy  is  slight,  and  if  S.  took  substantially  less  execu¬ 
tion  time  than  S.,  then  B.  is  to  be  preferred  to  B.,  for  it  provides  greater 
speed  with  comparable  accuracy. 


In  order  to  compare  S.  with  S.,  locate  m  points  of  time,  t , ,  t^ ,  .  .  .  ,  t  ,  at 

which  solutions  appear  in  both  sets  of  output.  Assuming  B^  more  restrictive 

than  B.,  a  good  measure  of  the  discrepancy  at  t  =  t,  is  given  by: 

J  ^ 

15 

-  -  ^ 

15 


D, 


n  =  1 


in 


where  s.  and  s.  are  the  n  functions  in  S.  and  S.  respectively.  If  D,  is 
in  jn  1  J  ^  ^  k 

small  for  k  =  1,  .  .  .  ,  m,  then  B.  and  B.  provide  about  the  same  accuracy. 

1  J 

Be  sure,  when  seeking  to  optimize  GLB  and  LUB,  that  the  initial  value  of 

the  step  size,  h^,  is  reasonable.  If  h^  is  too  big,  the  accuracy  will  be  poor 

no  matter  how  small  an  upper  bound  is  used.  An  initial  step  size  of  0,0  5  sec 

is  usually  satisfactory  for  the  WCHA  program.  However,  as  a  precaution,  it 

-7  -  5 

might  be  wise  to  begin  with  h^  =  0.  001  sec,  GLB  =  10  ,  LUB  =  10  ,  and  see 

how  the  step  size  behaves  over  an  integration  period  of,  say,  60  sec.  Having 

thus  determined  what  range  of  step  sizes  is  appropriate  for  the  system  at  hand, 

select  h  from  the  lower  end  of  that  range, 
o  ^ 


33 


Summary 


The  Adams- Moulton  subroutine,  AM(H,  T,  Y,  FY,  GLB,  LUB),  computes 
the  integral  curves  for  the  system  X. 

A.  Arguments 

H,  REAL  8;  input;  initial  step  size. 

T,  REAL  8;  input;  T  =  t^^ 

Y(15),  REAL  8;  input;  Y  = 

FY(15,  5),  REAL  >1^  8;  input;  FY(J,  K)  =  fj(t.  _  J  =  1,  15, 

K  =  1,3.  For  K  =  4,  5,  FY  is  computed  within  the  subroutine. 
GLB;  input;  greatest  lower  bound  on  error, 

LUB,  REAL;  input;  least  upper  bound  on  error. 

B.  Common 

/CHID0T/YD0T(15);/C0EF/P(4),  C(4) 

REAL  >:<  8:  YDOT,  P,  C 

C.  Subroutines  Called 

DERIV;  evaluates  the  derivatives. 

OUTPUT;  writes  out  the  solutions,  and  determines  when  to  termi¬ 
nates  the  simulation. 

RK;  restarts  AM  when  step  size  changes. 

Auxiliary  Subroutines 

I.  ANGMOM 

Subroutine  ANGMOM  (Y,  EL)  computes  I^,  the  angular  momentum  of  the 
satellite.  L  is  sometimes  needed  for  the  calculation  of  n  in  subroutine  TORQUE. 
In  the  special  case  of  torque-free  motion  (ji  -  X  “  angular  momen¬ 

tum  should  be  perfectly  constant;  the  accuracy  of  the  integration  can  then  be 
checked  by  writing  out  _L  and  observing  its  variance. 

A.  Arguments 

Y(15),  REAL  >1'  8;  input;  Y  = 

EL(3,  3),  real  8;  output;  EL(i,j)  =  L^^^,  i  =  1,  3,  j  =  1,  3. 
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B.  Common 

/MTM/OMEGA  (3,3),  IINV  (3,3),  DETI 
REAL  *  8:  OMEGA,  IINV,  DETI 

II.  DERIV 

Subroutine  DERIV  (Y,  YDOT)  evaluates  the  derivative  functions. 

A.  Arguments 

Y(15),  REAL  '•'  8;  input;  when  the  calling  program  is  AM,  either 

Y  =  ^  ^  “  Yi  4.  1  ^  4  when  the  calling 

program  is  RK,  Y  is  either  or  one  of  the  Runge-Kutta 
variables  ZJ_(t^),  Z^(t^),  or  Z3  (t.) 

YDOT(15),  REAL  8;  output;  YDOT  =  dy^/dt,  where  y_  input 

vector. 

B.  Common 

/INIT/EN(3);  /L0GIC/CH(3),  TRQ,  WHL;  /PARAM/IA(3),  R(3), 
KAY(3),  W,  W3. 

LOGICAL:  CH,  TRQ,  WHL 

REAL  8;  EN,  LA,  R,  KAY,  W,  W3. 

C.  Subroutines  Called 

RCTN;  evaluates  the  torque  of  the  reaction  wheel  on  the  satellite. 
RCTN  is  called  only  if  WHL  =  .  TRUE.  . 

TORQUE;  evaluates  the  external  torque  on  the  satellite.  TORQUE 
is  called  only  if  TRQ  =  .  TRUE.  . 

III.  OUTPUT 

Subroutine  OUTPUT(T,  Y,  HLV,  DEL,  NSTEP,  TERM)  is  a  user-  supplied 
output  routine,  called  by  AM  after  every  integration  step. 

A.  Arguments 

T,  REAL  8;  input;  T  =  t^ 

Y(15),  REAL  *  8;  input;  Y=x(ti)- 
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HLV,  INTEGER;  input;  number  of  times  h  has  been  halved. 

DEL,  INTEGER;  input;  number  of  times  h  has  been  doubled. 

NSTEP;  input;  number  of  Adams- Moulton  integration 
steps  which  have  been  taken. 

T  ERM,  LOGICAL;  output;  setting  TERM  to  .TRUE,  will  terminate 
the  simulation. 

B.  Common 

/OTPT/WR,  EPS,  TBOUND 

REAL  *  8:  WR 

REAL  *  4;  EPS 

/OTPT/  appears  on  WCHA  for  use  in  OUTPUT.  TBOUND  is  an 
upper  limit  on  T;  i.  e.  (T.  GE.  TBOUND)  can  be  used  as  a 
termination  criterion.  WR  and  EPS  may  be  used  to  obtain 
output  every  WR  sec  of  simulation  time.  The  following  in¬ 
structions  will  cause  the  solutions  to  be  written  out  if  and 
only  if  |T-n  *  WR  |  ^  EPS,  n  =  0,  1,  2,  .  .  . 

TMWR  =  DMOD  (T,  WR) 

IF  (TMWR.  GT.  EPS.  AND.  (WR-TMWR).  GT.  EPS)  GO  TO  10 

WRITE  (6,  20)  T,  Y 
20  FORMAT  (  •  •  •  ) 

10  GONTINUE 

Select  an  EPS  small  enough  to  constrain  the  output,  yet  not  so  small  as  to 
prohibit  it. 

In  choosing  EPS,  it  is  helpful  to  know  the  range  and  average  size  of  h. 
To  find  the  average  size,  write  out  NSTEP,  HLV,  and  DBL  at  the  end  of 
simulation,  and  divide  the  period  of  integration  by  the  total  number,  N,  of 
steps  taken;  N  is  given  by  NSTEP  +  3'!'(HLV  +  DBL  +  1).  To  get  an  idea  of 
the  range  of  h,  make  a  trial  run  with  a  small  TBOUND,  a  small  WR,  and  a 
large  EPS.  Another  way  to  study  the  behavior  of  h  is  to  write  out  HLV  and 
DBL  periodically  during  the  integration. 
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IV.  MXPROD 

Subroutine  MXPROD  (L,  M,  N,  MXL,  MXR,  MXP)  calculates  matrix  prod¬ 
ucts.  MXP  =  {MXL)(MXR),  where  MXL  is  L  by  M,  MXR  is  M  by  N,  and 
MXP  is  L  by  N.  All  the  matrices  are  REAL  8  arrays. 

V.  RCTN 

Subroutine  RCTN  (Y,  REAC)  is  a  user- supplied  routine  which  calculates 
the  torque  of  the  reaction  wheel  on  the  satellite.  The  form  of  the  routine  de¬ 
pends  on  the  control  law  governing  the  reaction  wheel. 

A.  Arguments 

Y(15),  REAL  8;  input;  Y  =  Y(t.) 

REAC,  REAL  8;  output;  reaction  torque  on  the  satellite 

VI.  TORQUE 

Subroutine  TORQUE  (Y,  EN)  is  a  user- supplied  routine  which  calculates 
n,  the  external  torque  in  body-fixed  coordinates.  The  form  of  the  routine  is 
determined  by  the  control  law  governing  the  external  torque. 

A.  Arguments 

Y(15),  REAL  8;  input;  Y  =  Y(t.). 

EN(3),REAL  8;  output;  EN  =  n(t.). 

Main  Program 

In  the  main  program,  WCHA,  all  the  input  data  is  read  in,  and  all  the 
COMMON  variables  are  initialized.  The  input  loop  has  the  form: 

1  READ(5,  DATA) 

IF  (H0.  EQ.  0.  )  GO  TO  6 
READ(5,  7)  TITLE 

7  FORMAT  (20A4) 

[  Integrate.] 

GO  TO  I 

6  CONTINUE 
STOP 
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A.  namelist/data/ 


All  the  variables  (except  TITLE)  are  read  in  on  a  single  list,  /DATA/. 
The  NAMELIST  variables  are; 

REAL  8;  initial  step  size.  Last  data  card  should  set  to  zero. 

Tj2^,  REAL  8;  initial  value  of  T,  not  destroyed. 

Y^{15),  REAL  8;  initial  value  of  Y,  not  destroyed. 

GLB;  greatest  lower  bound  on  error. 

LUB,  REAL;  least  upper  bound  on  error, 

WR,  REAL  8,  and  EPS(REAL  4);  output  can  be  obtained  at  desired 
intervals  by  writing  out  only  when  |T  -  n  *  WR  |  ^  EPS,  for 

n  =  0,  1,  2, . ,  , 

TBOUND;  upper  bound  on  T.  TBOUND-  T0  =  period  of  integration. 
TRQ,  LOGICAL;  TRQ  =  .  TRUE,  if  there  is  a  non- zero  external  torque, 
n,  on  the  satellite.  If  set  TRQ  to  .  FALSE.  . 

IA(3),  REAL  8;  IA(j)  =  L,  the  moment  of  inertia  of  the  satellite  about 
its  principal  axis. 

R(3),  REAL  8;  if  there  is  a  spherical  damping  rotor  on  the  axis, 

R(j)  =  R.,  its  moment  of  inertia.  If  there  is  no  rotor  on  the 
th 

j  axis,  set  R(j)  to  zero. 

KAY(3),  REAL  8;  if  there  is  a  spherical  damping  rotor  on  the 

axis,  KAY(j)  =  k.,  its  damping  constant.  If  there  is  no  rotor 
th  . 

on  the  j  axis,  set  KAY(j)  to  zero. 

W,  W3,  REAL  8;  if  the  satellite  has  a  reaction  wheel,  with  moments 
of  inertia  W^,  W^,  W^,  then  W  =  =  W^,  and  W3  =  W^.  If 

there  is  no  reaction  wheel,  set  W  and  W3  to  zero. 

B.  TITLE 

A  title  of  up  to  128  characters  is  printed  at  the  head  of  each  set  of  out¬ 
put.  It  is  read  in  as  an  alphameric  array  TITLE  (32),  Since  one  data  card  con¬ 
tains  only  80  characters,  two  cards  must  always  be  supplied  for  the  title.  If 
only  one  card  is  supplied,  TITLE(21),  TITLE(22),...,  TITLE(32)  will  be  read 
from  the  following  NAMELIST  card,  with  dire  consequences. 
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C.  COMMON 


The  COMMON  blocks  are  listed  below  with  the  subroutines  in  which 
they  appear. 

/CHIDOT/YDOT  (15);  RK,  AM 
/COEF/P  (4),  C(4);  AM 
/INIT/EN(3);  DERIV 
/LOGIC/CH  (3),  TRQ,  WHL;  DERIV 
/MTM/OMEGA  (3,  3),  IINV  (3,  3),  DETI;  ANGMOM 
/OTPT/WR,  EPS,  TBOUND;  OUTPUT 
/PARAM/IA  (3),  R  (3),  KAY  (3),  W,  W3;  DERIV 
LOGICAL;  CH,  TRQ,  WHL 

REAL  8:  YDOT,  P,  C,  EN,  OMEGA,  IINV,  DETI,  WR,  lA,  R,  KAY,  W,  W3 
Summary 

To  use  the  attitude  stabilization  program,  read  in  the  input  data  as  de¬ 
scribed  above,  and  supply  subroutines  TORQUE,  WHEEL,  and  OUTPUT.  If 
TRQ  =  .  FALSE.  ,  TORQUE  will  not  be  called,  and  if  W  =  0,  WHEEL  will  not 
be  called. 

Before  making  any  modifications  in  the  program,  study  the  interdepend¬ 
ence  of  the  subroutines  involved,  and  double-check  the  specifications  of  the 
variables.  The  pertinent  information  is  tabulated  in  the  above  subroutine 
profiles. 
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